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We present a state-of-the-art primordial recombination code, HyRec, including all the physical 
effects that have been shown to significantly affect recombination. The computation of helium 
recombination includes simple analytic treatments of hydrogen continuum opacity in the He I 2 1 P° — 
1 1 5 line, the He I] 2 3 P° — US line, and treats feedback between these lines within the on-the-spot 
approximation. Hydrogen recombination is computed using the effective multilevel atom method, 
virtually accounting for an infinite number of excited states. We account for two-photon transitions 
from 2s and higher levels as well as frequency diffusion in Lyman-a with a full radiative transfer 
calculation. We present a new method to evolve the radiation field simultaneously with the level 
populations and the free electron fraction. These computations are sped up by taking advantage of 
the particular sparseness pattern of the equations describing the radiative transfer. The computation 
time for a full recombination history is ~ 2 seconds. This makes our code well suited for inclusion 
in Monte Carlo Markov chains for cosmological parameter estimation from upcoming high-precision 
cosmic microwave background anisotropy measurements. 



I. INTRODUCTION 

Until recently, primordial recombination was consid- 
ered one of the few solved problems in astrophysics and 
cosmology. The seminal works of Peebles [JJ and Zel- 
dovich et al. [2] in the 1960s established that hydrogen 
recombination did not proceed in Saha equilibrium, and 
that two-photon decays from 2s are critical to the re- 
combination dynamics because of the very low escape 
rate of Lyman-a photons. They provided a simple ef- 
fective three-level atom model to compute primordial 
hydrogen recombination histories. With the advent of 
high-precision cosmic microwave background (CMB) ex- 
periments such as WMAP [3J and Planck [2], it has be- 
come clear that these early calculations are not suffi- 
cently accurate for an unbiased estimate of cosmologi- 
cal parameters [5H7] . Uncertainties in the recombination 
history indeed propagate to the visibility function and 
ultimately to the predicted CMB temperature and po- 
larization anisotropies. 

These considerations have motivated Seager et al. 
[U [9] to extend Peebles' effective three-level atom model 
to a multilevel atom calculation. They showed that ac- 
counting for excited states of hydrogen leads to a speed 
up of recombination at late times. Their commonly used 
recombination code RecFast approximately reproduces 
these results by solving an effective three-level atom 
model with an artificially enhanced recombination co- 
efficient. While this model is sufficiently accurate for 
current CMB data analysis, it does not meet the ~ 0.1% 
accuracy target for Planck. 

Several physical effects have since then been shown 
to significantly affect hydrogen and helium recombina- 
tion. First, the multilevel computations of Seager et 
al. assumed statistical equilibrium among the angular 
momentum substates of a given energy shell of hydro- 
gen. At late times, this assumption breaks down, and 
an accurate multilevel atom computation should resolve 



the angular momentum substates [TO] ITT] . Whereas it 
is a straightforward conceptual generalization of previ- 
ous works, this problem can be computationally chal- 
lenging. Several works have tackled the "high-n prob- 
lem" , including increasingly larger numbers of excited 
states of hydrogen to reach sufficient accuracy [TfJlfT3] . 
The "standard" multilevel atom (MLA) method used in 
these works requires solving for the abundances of all the 
excited states at every timestep, which makes the compu- 
tation very time-consuming. In a recent paper (Ref. |14j , 
hereafter "Paper I"), we have introduced a new method 
of solution for the multilevel atom problem, which allows 
to only solve for the populations of a few excited states 
(typically, 2s and the low-lying p states), provided one 
uses precomputed effective bound-bound and bound-free 
transition rates, which contain all the information about 
the highly excited states of hydrogen. This method al- 
leviates the computational difficulty associated with the 
highly excited states of hydrogen and is key to the speed 
of the recombination code presented here. 

Another important aspect of the recombination prob- 
lem is that of radiative transfer in the vicinity of the 
Lyman-a line. In its early stages, hydrogen recombina- 
tion is mostly controlled by the slow escape (via redshift- 
ing) of photons from the Lyman-a line and the rate of 
two-photon decays from the 2s state. Accurate values for 
these rates require treatements of the radiation field that 
go beyond the simple Sobolev approximation |8, 15 . Im- 
portant corrections include feedback from higher-order 
lines [TH] [T7], time-dependent effects in Ly-a [THj, and 
frequency diffusion due to resonant scattering [T9Tl21| . 
An accurate 2s — Is two-photon decay rate also requires 
following the radiation field to account for stimulated de- 
cays [22] and absorption of non-thermal photons [22] [21] . 
Dubrovich & Grachev [IS] suggested that two-photon 
transitions from higher levels may have a significant ef- 
fect on the recombination history. Later computations 
confirmed this idea [26] , and provided an accurate treate- 
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ment of radiative transfer in the presence of two-photon 
transitions, as well as a solution for the double-counting 
problem (which arises for resonant two-photon transi- 
tions, already included in the one-photon treatment as 
"1+1" transitions) 1231127] . 

The accuracy requirement is less stringent for pri- 
mordial helium recombination, as it is completed by 
z ~ 1700, much earlier than the peak of the visibility 
function. Corrections at the percent level are still impor- 
tant, and several works have been devoted to the problem 
5, 8, 25, 28 35j. The most important effect is continuum 
opacity in the He I 2 1 P° — l 1 ^ line due to photoioniza- 
tion of neutral hydrogen, which requires a detailed line 
+ continuum radiative transfer analysis [30l IS , 34 . The 
inclusion of the intercombination line He I] 2 3 P° — l 1 ^ 
is also significant. 

Several other processes have been investigated and 
shown not to be significant for CMB anisotropies, for ex- 
ample the effects of the isotopes D and 3 He [TH [TTJ [33J , 
lithium recombination |36j . quadrupole transitions |12j . 
high-order Lyman line overlap |15j . and Thomson scat- 
tering JXSJ EI] • Collisional processes are negligible for he- 
lium recombination |33] ; for hydrogen recombination, col- 
lisional corrections appear to be small |37j . but whether 
they are truly negligible is still under investigation. 

Previous works have all concentrated on one or a few 
aspects of the primordial recombination problem. Pro- 
ducing a complete and fast recombination code has so 
far been hindered by the computational burden previ- 
ously associated with the high-n problem. Given that 
this problem is now solved, and that it seems that the 
main radiative transfer effects have now all been identi- 
fied, it is timely to deliver a single code that computes 
an accurate hydrogen and helium recombination history 
and incorporates all the relevant physics. The purpose of 
this paper is to introduce our new recombination code, 
HyRec, which is publicly available^ and can compute a 
highly accurate recombination history (with errors at the 
level of a few times 10 for helium recombination and a 
few times 10 -4 for hydrogen recombination) in only ^2 
seconds on a standard laptop. Our code does not account 
for collisional transitions in hydrogen, as their rates are 
poorly known. When accurate rates are available and if 
collisional transitions are shown to significantly impact 
recombination, we will update our code with the appro- 
priate effective rates. 

Recently, a similar work has been carried out by 
Chluba & Thomas [35J , also relying on the effective ML A 
method presented in Paper I. The code they present in- 
cludes the same physics as ours. The main difference is 
the treatment of radiative transfer. In Ref. [35], an "or- 
der zero" recombination history is first computed, with 
a simple treatment of radiative transfer. The radiative 
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transfer equation is then solved, given this order zero his- 
tory. Corrections to the net decay rates to the ground 
state are then evaluated, and used to compute a corrected 
recombination history. This procedure can in principle 
be iterated, but because the corrections are small, it is 
essentially converged after one iteration. Our solution, 
on the other hand, is non-perturbative, in the sense that 
we solve simultaneously for the radiation field and the 
recombination history. A detailed code comparison is in 
progress, and a full error budget will be presented once 
it is completed. 

This paper is organized as follows. In Sec. [H] we re- 
view the effective three-level atom model and discuss its 
limitations. We then review the standard MLA compu- 
tation and describe the effective MLA method in Sec. Mil 
We show that weak transitions to the ground state from 
excited states with n > 3 can in fact be accounted for 
almost exactly with an effective four-level atom model. 
Two-photon processes and frequency diffusion are for- 
mally described in Sec. IV In Sec. [V] we present our 



numerical solution for the radiative transfer equation. 
We use a new method of solution, extending that of 
Ref. [UJ to account for frequency diffusion, that allows 
to solve for the atomic populations and the radiation 
field simultaneously. We describe our treatment of he- 
lium recombination in Sec. IVII We conclude in Sec. IVHI 
Appendix [A] demonstrates some relations satified by the 
effective rates, Appendix [B] describes how we extrapo- 
late the effective rates to an infinite number of excited 
states, Appendix [C] describes our ordinary differential 
equation (ODE) integrator, and Appendix |d] derives an 
analytic expression for the post-Saha expansion used at 
early times in hydrogen recombination. 

Throughout this paper we use a flat background 
ACDM cosmology with T = 2.728 K, n b h 2 = 0.022, 
fl m h 2 = 0.13, n A h 2 = 0.343, y Ho = 0.24 and N^ c g = 
3.04. 



II. HYDROGEN RECOMBINATION: 
OVERVIEW 

A. The effective three-level atom model 

The basic process of primordial hydrogen recombina- 
tion was already well understood in the late sixties. The 
seminal papers by Peebles [T] and Zeldovich et al. E] 
established the following picture. Direct recombinations 
to the ground state are highly inefficient, as they pro- 
duce photons that can immediately ionize another hy- 
drogen atom. Electrons and protons can therefore re- 
combine efficiently only to the excited states of hydrogen. 
This situation is familiar in the study of the interstellar 
medium: it is referred to as "case-B" recombination (see 
e.g. Ref. [39 ). Once they have recombined to one of 
the excited states of hydrogen, electrons "cascade down" 
to the n — 2 shell, on a much shorter timescale than 
the overal recombination timescale. Denoting nn the to- 
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tal number density of hydrogen, x e — n e /n-a the free 
electron fraction and x% = i-H(n=2)/ n H the fraction of 
hydrogen in the excited state, the effective rate of re- 
combinations to n — 2 shell can be written: 



3'2 



-x e = n H x 2 e a B {T m ) - x 2 f3 B (T T ), (1) 



where T m is the matter temperature, locked to the ra- 
diation temperature T r by Thomson scattering at most 
times during recombination, a B is the case-B recombina- 
tion coefficient, and f3 B is the corresponding photoioniza- 
tion rate, which can be obtained from a B by the principle 
of detailed balance: 



0B(T r ) 



n H a B (T m = T r ), 



where we have defined 
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where [i e is the reduced mass of the electron-proton sys- 
tem. 

Once they have reached the n = 2 shell, electrons can 
reach the ground state by emitting a Lyman-a photon 
from the 2p state. Due to the very high optical depth of 
the Lyman-a transition, emitted photons will however al- 
most certainly be reabsorbed by another atom. The way 
out of this bottleneck is for photons to redshift below 
the Lya resonant frequency due to cosmological expan- 
sion. The net rate of decays to the ground state from the 
2p state is then just the rate at which photons redshift 
across the line and escape reabsorption: 



x u 



2p 



~ X 2p\ls 



i? LyQ (x 2p - 3x ls e- £21 / T ') , (4) 



where xu is the fraction of hydrogen in the ground state 
and the second term accounts for Lya absorptions and is 
obtained by detailed balance. The rate of escape of Lya 
photons is given by: 
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Eqs. |4jj5]) can be derived in the Sobolev approximation, 
in the limit of large Sobolev optical depth (see for exam- 
ple Ref. QU). 

The escape rate of Lya photons is comparable to the 
rate of the slow two-photon decays from the 2s state, 
A2s,is ~ 8.22 s _1 , and the latter process must therefore 
be accounted for. The net rate of two-photon decays from 
the 2s state is: 
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xi s e 



-E 2 i/T r 



(6) 



where the second term accounts for two-photon absorp- 
tions and can be obtained by a detailed balance argu- 
ment. Due to the strong thermal radiation bath, the 
excited states of hydrogen are near Boltzmann equilib- 
rium with each other, x 2p = 3x 2s = (3/4)x2- The rate of 



change of the population of the n = 2 shell due to decays 
to the ground state is therefore: 
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The last step is to realize that the atomic rates, even for 
the slow 2s — > Is decays or the slow escape out of the Lya 
resonance, are many orders of magnitude larger than the 
overall recombination rate, which is of the order of (10 
times) the Hubble expansion rate, that is ~ 10~ 13 — 1CP 12 
s" 1 . The population of the n = 2 shell can therefore be 
obtained to high accuracy in the steady-state approxima- 
tion, i.e. assuming that the rate of recombinations to the 
n = 2 shell equals the rate of transitions to the ground 
state: 



x 2 = x 2 
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We can therefore solve for x 2 and obtain: 



x 2 



riB\xla B + (3-R L ya + A 23 ,is) xi s e 

+ jRhya + |A 2s! ls 



-E 2 i/T t 



(8) 



(9) 



From Eq. wc then obtain the rate of change of the 
free electron fraction: 



x e = -C (riB\xla B - Axi s /3 B e 
where the Peebles C-factor is given by 
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(10) 



(11) 



As noted by Peebles, this factor represents the proba- 
bility that an atom initially the n = 2 shell reaches 
the ground state before being photoionized. Note that 
we could have obtained the same equation starting from 
x e = —xu = —(xu\2p + iisbs) (this is because we have 
set x 2 = 0). 

At all relevant times during the epoch of hydrogen re- 
combination, x 2 -C 1, and therefore x\ s = 1 — x e . If 
matter and radiation temperatures are set to be equal, 
Eq. ( 10 1 is therefore a simple ordinary differential equa- 



tion for x e , that can be easily integrated. A simple im- 
provement is to also explicitly follow the matter temper- 
ature evolution, which is determined by the Compton 
evolution equation: 



-2HT„ 



&(TTa r Tfx e (T r - T m ) 
3(1 + / He + x e )m e c 



(12) 



where <7t is the Thomson cross-section, a r is the radiation 
constant, m e is the electron mass and /n e is the He:H 
ratio by number of nuclei. 

The simple yet insightful picture presented here is 
known as the effective three-level atom model. It pro- 
vides a good approximation for the recombination prob- 
lem. It is however not sufficiently accurate for high- 
precision cosmology. 
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FIG. 1. Peebles C-factor [Eq. ( 11 1] and ratio of the population 
of the n = 2 shell to its value in Saha equilibrium with the 
continuum, as a function of redshift. 



B. Hydrogen recombination phenomenology 

We show in Fig. [T] the evolution of the Peebles C-factor 
and the population of the n — 2 shell relative to its value 
in Saha equilibrium with the continuum, as a function of 
redshift, for a standard recombination history. We can 
see that there are two distinct regimes. 

At early times (z > 1000), electrons in the n = 2 
shell have a high probability of being photoionized, and 
the C-factor is much smaller than unity, C < 1. As a 
consequence, the population of the n = 2 shell is very 
close to Saha equilibrium with the continuum, 



x 2 ~ x 2 



Saha 



9e 



(13) 



The rate of change of the free electron fraction is then 
approximately equal to the rate of decays from the n = 2 
shell: 



x e {z> 1000) wi 2 1 (x 2 =x 2 \ Saha ) 



(14) 



During that period, the recombination rate is therefore 
virtually independent of the exact value of the recombi- 
nation coefficient, but is strongly dependent on the small 
net decay rate from the n = 2 shell to the ground state. 
This is usually referred to as the "n = 2 bottleneck" and 
has motivated abundant work on radiative transfer in the 
vicinity of the Lyman transitions [T5rH5I l2"0H2"4l [2^1127] , 
The result from this series of papers is that to the level 
of accuracy required by Planck, Lyman transitions up to 
Ly7 must be included, properly accounting for feedback 
between them. In addition, the radiation field must be 
solved for with a radiative transfer calculation, account- 
ing for two-photon transitions and frequency diffusion in 
the Lyman-a line. We consider all these effects in Sec- 
tions |IV] and El 

At late times (z < 700), C s» 1, and the n = 2 shell is 
no longer in Saha equilibrium with the continuum (note 
that it is not in Boltzmann equilibrium with the ground 
state either, as the rate of recombinations to the n = 2 
shell dominates over the net rate of two-photon or Ly-a 



absorptions from the ground state). The free electron 
fraction is many orders of magnitude above the value 
it would have in Saha equilibrium because of the slow 
recombination rate. In that case, the second term in 



Eq. ( 10 ) is negligible and the evolution of the free electron 



fraction becomes: 



x e (z < 700) 



(15) 



As we can see, the evolution of the free electron frac- 
tion is then virtually independent of the rate of decays 
to the ground state from the n — 2 shell, but is highly 
sensitive to the exact value of the effective recombina- 
tion coefficient. Moreover, the assumption that the ex- 
cited states are in Boltzmann equilibrium with each other 
brakes down at late times because of the decrease in the 
radiation temperature. An accurate recombination rate 
at late times can only be obtained in a full multilevel 
atom calculation, that accounts for (possibly stimulated) 
bound-bound and bound-free transitions between all - 
at least, a large number of - the excited states of hy- 
drogen [SI [TTI - H5] . We will review the multi-level atom 
calculations in Sec. IIII1 

Of course, at intermediate redshits 700 < z < 1000 
both the exact recombination rate and the rate of decays 
to the ground state are important and should be carefully 
accounted for. 



III. THE MULTI-LEVEL ATOM 

In this section we first present the "standard" multi- 
level atom (MLA) method [Hj, then review the effective 
MLA (hereafter EMLA) method of solution that we pre- 
sented in Paper I and that allows for a fast computation 
of recombination histories. 



A. The standard multi-level atom method 

The effective three-level atom equations may be easily 
generalized to account for an arbitrarily large number of 
excited states of hydrogen (in practice, one must of course 
impose a cutoff). We denote x n i the fractional abun- 
dance of hydrogen atoms in the excited state with prin- 
cipal quantum number n and angular momentum quan- 
tum number I. The generalization of Eq. ([!]) is then, for 
n > 2: 



nnx 2 e a nl (T m , T T ) - x n i/3 n i(T,), (16) 



where a n i(T m , T r ) is the recombination coefficient to the 
excited state nl, including stimulated recombinations, 
and /3 n i(T T ) is rate of photoionizations from nl by black- 
body photons. 

The effective three-level atom model does not ac- 
count for bound-bound transitions between excited states 
(except for instantaneous spontaneous decays that ulti- 
mately lead to the n = 2 shell). Transitions between the 
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nl and n'l' states (with n, n' > 2) change their popula- 
tions at the rate: 



Xnl 



n'l' 



Xn'l' 



nl 



X n 'l'Rn'l',nl(T T ) — X „/ R n l,n' I' (T T ) , (17) 



where the bound-bound transition rate from nl to n'l 1 , 
Rni,n'l'(Ti), is the rate of absorptions of blackbody pho- 
tons resonant with the transition if n < n' and the rate 
of spontaneous and stimulated decays if n > n'. We give 
explicit expressions of the bound-bound and bound free 
rates and explain how we compute them in Paper I. 

Finally, the rate of decays to the ground state from the 
2s state are given by Eq. Q (we will see how to make 



this rate more accurate in Sec. IV A). In the Sobolev 
approximation, the net decay rate from the np states to 
the ground state is given by a generalization of Eq. Q , 
accounting for feedback between optically thick Lyman 
lines: 



xi s 



np 



where Rh yn = {^Lya/ ALyn) 3 Rhya is the rate at which 
photons rcdshift out of the Lyman-n line (with the con- 
vention that Ly-2 is Lya), and /+, is the photon occu- 
pation number incoming on the blue side of the Ly-n 
transition. If no radiative processes affect the radiation 
field between neighboring Lyman lines, then 



fnp( z ) ~ fn+l,p( z )) 

where the earlier redshift z' is given by 

/ _ ^Lyn 



A 



Ly(n+1) 



(l + z)-l. 



(19) 



(20) 



In the optically thick limit which is valid here, the photon 
occupation number redward of the Ly-n line is given by 

I np = x np/{3xis)- 

For excited states nl other than 2s and np, not radia- 
tively connected to the ground state, we have i n jL = 0. 

The recombination history can then be computed by 
evolving simultaneously the system of differential equa- 
tions: 



Xnl — X n l 



E 



Xnl 



n ' I ' 



Xnl 



Is 



n'>2,l' 
-Xu = X 2 s\ ls " 



(21) 
(22) 



where in the last equation we used x e = 1 — xi s , valid 
as the fractional abundance of hydrogen in the excited 
states is always much less than unity. The atomic transi- 
tion rates are many orders of magnitude larger than the 
overall recombination rate (of the order of the Hubble 
rate). A highly accurate approximation therefore consists 
in first solving for the populations of the excited states in 
the steady-state approximation. This first step amounts 
to solving a large system of linear algebraic equations. 



The populations X2 S , x np can then be used in Eq. ( 22 1 to 
evolve the free electron fraction. 

This generalization of the three-level atom model is 
relatively straightforward conceptually. Its practical im- 
plementation is, however, very time-consuming. It re- 
quires solving a very large system of algebraic equations 
at each time step (or evolving the same number of stiff 
differential equations) . If one accounts for excited states 
up to principal quantum number n max , then the num- 
ber of equations is N = n max (n max + l)/2, which, for 
n max > 100, exceeds several thousands. Furthermore, 
the computational cost of an exact linear system solu- 
tion scales as 0(N 3 ), although this can be significantly 
sped up by using the sparseness of the system due to se- 
lection rules [T2] and iterative solution techniques [37] . In 
the following section, we review the much more efficient 
yet exactly equivalent effective MLA method. 



B. The effective multi- level atom method 

1. General description 

The effective multilevel atom method, described in Pa- 
per I, relies on three aspects of the primordial recombi- 
nation problem. First, the timescales for transitions out 
of the excited states are much shorter than the overall re- 
combination timescale - this property is used when solv- 
ing for the populations of the excited states in the steady- 
state approximation in the standard MLA method. This 
allows us to factor all the nearly instantaneous transi- 
tions involving the "interior" excited states (which are 
not radiatively connected to the ground state) into effec- 
tive transitions into and out of the smaller set of "in- 
terface" states which are radiatively connected to the 
ground state. Secondly, all bound-bound and bound- 
free transitions for which the lower state is an excited 
state are optically thin, and therefore do not distort the 
ambient blackbody radiation field in the vicinity of the 
corresponding frequencies. All the transition rates be- 
tween "interior" states therefore only depend on the ra- 
diation temperature T r (as well as atomic physics con- 
stants). This translates into a simple dependence for the 
effective rates, which are functions of the matter (through 
recombinations of thermal electrons and protons) and ra- 
diation temperatures only. If collisional transitions are 
included, they will depend additionally on the free elec- 
tron (of equivalently the free proton) abundance. Finally, 
the set of "interface" states that need to be considered 
is small. In principle, 2s and all the p states need to 
be considered as "interface" states; however, in practice 
only the lowest order Lyman transitions significantly af- 
fect the recombination rate, explicitly only Lya, Ly/3 and 
Ly7 p31[IH]. In addition, when considering two-photon 
transitions from higher levels, one should in principle add 
the ns and nd states as "interface" states. However, only 
two-photon transitions from 2s, 3s and 3c?, are important 
(and 4s, Ad at the level of a few 10" 4 ) [H [27J. This 
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means that one needs to pretabulate only a few func- 
tions of temperature that fully account for the multilevel 
structure of hydrogen. These tabulated effective rates 
can then be interpolated when computing a recombina- 
tion history with an effective few-level atom model. As 
we will show in Sec. |III B 2[ we can in fact further sim- 
plify the problem to an effective four-level atom model 
with virtually no loss of accuracy. 

Following Paper I, we denote Ai(T m ,T r ) the effective 
recombination coefficient to the "interface" state i, Bi(T r ) 
the effective photoionization rate from this state, and 
lZij(T r ) the transfer rate from the interface state i to 
the interface state j (the dependences are valid in the 
purely radiative case; when collisions are included, all 
effective rates depend on T m ,T r ,n e ). They are obtained 
as follows: for the effective recombination coefficients, 



(23) 



A' 



for the effective ionization coefficient 
Bi = A + £ R ^ kP ' 



Ki 



K 



and for the effective inter-state transition rates, 



(25) 



K 



Here K is a general index for "interior" states, P\^ is 
the probability that an electron initially in the "interior" 
state K ultimately reaches the "interface" state i, and 
Pjf is the probability that an atom initially in the state 
K ultimately gets photoionized. These probabilities are 
the solutions of the linear systems [14]: 
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where 
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T/< = ^ Rk.l + ^ Rk.i + Pk 



(26) 



(27) 



(28) 



is the inverse lifetime of the state K. 

The net decay rate from the interface state i to the 
ground state can always be expressed as a linear function 
of x, : 



I Is 



= XiR 



2,1S 



x\ s R\ t 



(29) 



where we emphasize with this notation that in general 
the net decay rates may depend in a complicated way on 
the current and past values of the free electron fraction, 
as well as on cosmological parameters - see for example 



Eqs. (18) and ^ for the net np — > Is decay rate. This 
contrasts with the bound-bound rates between excited 
states, which only depend on atomic constants and the 
radiation temperature. 

Once the effective bound-bound and bound-free rates 
for the interface states are tabulated, the recombination 
history can be computed by evolving the small set of 
ordinary differential equations: 



sRls 



R l , 



(30) 



and 



- (nuxlAi - Xl Bi) (31) 

i 

-iis = ^ (xisRis^i - XiRi^ , (32) 



(24) where Eqs. (31) and (32) are equivalent as the fractional 



abundance of excited hydrogen is very small. In practice, 
the population of the effective states can once again be 
solved in the steady-state approximation. Eq. ( 30 ) be- 



comes a system of a few algebraic linear equations, and 



in that case Eqs. ( 31 ) and ( 32 ) are mathematically equiv- 
alent since we set ii = [this can be seen by summing 



Eq. (30) over jL 



Equations (30) and (31 ), along with the definitions for 



the effective coefficients Eqs. (23-28), are strictly equiv- 



alent to the standard MLA equations presented in the 
previous section, as was derived in Paper I [Eq. (31) was 
not derived in that paper and we give a proof in Ap- 
pendix A 2 . The advantage of the new method is that 



the system of equations that need to be solved at each 
time step is much smaller, as only the low-lying s,p,d 
states need to be followed. In the following section we 
show that we can actually further reduce the problem to 
an effective four-level atom model. 



2. Further simplification: the effective four-level atom 



If we wish to follow "interface" states, then the sys- 



tem of equations (30), in the steady-state approximation, 
is a x system. Moreover, one needs to interpolate 
n* functions of 2 variables (the effective recombination 
coefficients - the effective photoionization rates are ob- 
tained by detailed balance), and n*(n* — l)/2 functions 
of 1 variable (half of the effective bound-bound rates, the 
other half being obtained by detailed balance) . Here we 
show how this system can be further reduced to a 2 x 2 
system involving only 2s and 2p, requiring only 2 func- 
tions of 2 variables and 2n* — 3 functions of one variable, 
with virtually no loss of accuracy. 

For now on we use the general index K for all states 
with principal quantum number n > 3. Even if some of 
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the states with n > 3 are radiatively connected to the 
ground state, one can still formally define the effective 
transition rates Eqs. (23), (24) and (25) for the 2s and 2p 



states because of the near-instantaneity of transitions out 
of the excited states. However, these coefficients do not 
have a simple temperature dependence anymore, and are 
therefore not well suited for fast interpolation. Indeed, 
the probabilities P % K and P%, where i = 2s, 2p, are still 
defined by Eqs. (26) and (27), but the inverse lifetime of 



the state K, Eq. (28), should now account for the net 



downward transition rate to the ground state, and one 
should make the replacement: 



K 



(33) 



where we define IV (T r ) is the inverse lifetime of the K-th 
interface state when transitions to the ground state are 
not included. In addition, we need to define additional 
effective transition rates with the ground state. For i = 
2s, 2p, we define: 



Hi 



R 



2,1s 



and 
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Is,: 



R 



IsA 



K 



/.Ris.kPk) 

K 



(34) 



(35) 



where the probabilities Pf? must satisfy the self- 
consistency relations: 



P, 



A 



r. 1 A 



R 



KAs 



r 



(36) 



A 



The standard MLA equations, in the steady-state ap- 
proximation for the excited states with n > 3, can then 
be shown to be exactly equivalent to the following set of 
equations: for the net rate of production of 2s, 

X 2s = x\n^Als + X2pR-2p,2s + X ls TZi St 2s 

-x 2s \B 2s + R-2 S ,2p + H2 S ,is) ; (37) 
for the net rate of production of 2p, 

X 2 p = xlnuA 2 p + X2sH 2s .2p + XlsR-ls^p 

— X 2 p {^2p + R-2p,2s + R-2p,ls^j ; (38) 

and for either the net recombination rate or the net rate 
of generation of H(ls), 



^ [x i B l - xlriYiAi} 

i=2s,2p 



i=2s,2p 



<£\s'R'\s,i ^{R%At 



(39) 
(40) 



The proof of equivalence is a simple generalization of that 
of Paper I and we do not reproduce it here. Note that 



we use the same notation for the effective rates indepen- 
dently of the number of "interface" states considered, but 
they obviously have a different meaning that should be 
clear from the context. 

The coefficients in the above system are in principle not 
simple functions of temperature anymore if one wishes 
to account for the transitions to the ground state from 
excited states with n > 3. We can nevertheless simplify 
their expressions with some minimal approximations. We 
start by noticing that for excited states nl with n > 3, 
the rate of spontaneous decays to the n',l ± 1 states 
with 1 < n' < n is much larger than the net decay 
rate to the ground state. For the 3p state for exam- 
ple, we find that in the Sobolev approximation for Ly/3 
decays, R 3p ,u/A 3pt2s < 8 x 1(T 6 for 200 < z < 1600, 
where A^p^s is the Einstein A-coefficient of the 3p —> 2s 
transition. Therefore, to an excellent accuracy (with rel- 
ative errors of order Rkas/^k), one can neglect Rkas in 
Eq. (33), and simply use IV (T r ) instead of IV wherever 
the latter appears. With this approximation, the rate co- 
efficients A2sj2pi $2s/2p an d R-2s,2 P are simply the usual 
effective rates computed in the case that only 2s and 2p 
are considered as interface states, and depend only on 
matter and radiation temperatures. We explain in Ap- 
pendix [B] how we obtain effective rates extrapolated to 

%iax — OO. 

We show in Appendix 



A3 



Eq. ( 36 ) , we can rewrite Eq. 



that using 
d34| as: 



T K (T r ) in 



Tli,ls = Ri,ls + E 



K 



i 9 —e- E 

Is e 

9% 



™'T*P l K (T r ), (41) 



where gK,gi are the statistical weights of the states K,i 
and Ek2 = Ek — E 2 is the energy difference between the 
state K and the n = 2 shell. 

In addition, the populations of the "weak" interface 
states Xk are sometimes required - for example the pho- 
ton occupation number depends on the populations of the 
s and d states, see Sect ion fVBj We show in Appendix |A4| 
that the following relation is verified for T m = T r : 



X 



A 



='A 



9K e -E K /T ip e K{Ti)x 2 e 



9e 

+ 9 -^e- E *"' T <P l K {T r )x t . (42) 

i=2s.2p 

In fact, T m < T r and the coefficient of x"l in the above 
equation should be slightly higher. In practice though, 
\T m /T T - 1| < 1% for z > 500, and for lower redshifts 
<C 1 and transitions to the ground state are unim- 
portant anyway, so the above equation is very accurate. 

We therefore only need to tabulate the additional 
2(rc» - 2) functions P| s (T r ) and P^- p (T r ) to account for 
— 2 "weak" interface states in addition to 2s and 2p 
(note that P£ = 1 - P|? - P^). 

In practice, we can further reduce the computational 
load by simply using P„ 2 f = P„ 2 ^ = P%> = 0, P n 2 ? = 



P 2p 
r nd 



P 2 p = 1. This amounts to assuming that for 



n > 3, ns and nd states are in Boltzmann equilibrium 




CUD 



FIG. 2. Schematic representation of the hydrogen atom, with 
the nomenclature used in this paper. Slow transitions from 
the "weak interface" states to the ground state are counted 
as transitions from the n = 2 state with which they are in 
equilibrium. 



with 2p, whereas np states are in Boltzmann equilibrium 
with 2s - see Eq. ( 42 ) , and rewriting transitions from 
n > 3 states as transitions from the n 



= 2 state to which 
they are tightly coupled. This is extremely accurate at all 
times when two-photon transitions significantly affect re- 
combination. The validity of this statement is somewhat 
weaker for the n — 4 states, which are somewhat in be- 
tween equilibrium with the 2s state and the 2p state (be- 
cause allowed transitions connect them to both states). 
However, as decays from 4p, 4s and Ad only marginally 
affect recombination anyway (at the level of a few 10 -4 ), 
and 2s and 2p are very close to equilibrium at the rele- 
vant times, this approximation is still very accurate. We 
explicitly checked that using the approximate valu es for 
the P % K instead of their exact values in Eqs. (35), (41) 
and ( 42 ) leads to maximum errors on the recombination 
history] Ax e | jx e < 3 x 10~ 5 . 

We illustrate the formulation adopted in this paper 
schematically in Fig. [2j 



The system of equations (37), (38) and (39) is just the 



extension of Peebles effective three-level atom to an effec- 
tive four-level atom, properly accounting for the non-zero 
radiation field, the nearly instantaneous multiple transi- 
tions between excited states, the fact that 2s and 2p are 
out of Boltzmann equilibrium, and possibly additional 
radiative transfer effects and decays from higher shells 
through the appropriate coefficients 1Z. Making the usual 
steady-state assumption for the excited states, we can 
first solve for x 2s and X2 P and then evolve x e . When us- 
ing the simple 2p O Is and 2s ■<-> Is transition rates of 
Sec. II A| the system is simple enough that we may write 



the function x e explicitly as an illustration: 

x e =- C 2s (n H x 2 e A 2s - x ls B 2s e~ E2l/T r 
- C 2p (n u x 2 e A 2p - 3x ls B 2p e~ E2l/Tr 



where the C-factors are given by 



C 



2 s 



2s 



K 



R 2 



2s^2p^T 2 



and 



C 



R 



Lya 



n 



2p 



A 2s 

2p^2s-p 2 



2p 



R, 



(43) 



(44) 



(45) 



2 P ^2s—p^ 



and where we have used the effective inverse lifetimes: 
r 2 s = B 2s + Tl 2s . 2p + A 2s ,i s and 

r 2p = B 2p + lZ 2p , 2s + Rhya- (46) 

In Fig. [3] we show the changes to the recombination his- 
tory resulting from an accurate effective multi-level com- 
putation, as compared to the effective three-level atom 
computation, using in both cases the simple decay rates 
to the ground state described in Sec II A For comparison, 
we also show the resulting changes when using an effec- 
tive three-level atom model with a fudge factor F = 1.14 
as in the code RecFast [§]. We checked that it is not 
possible to reproduce the correct effective MLA compu- 
tation with a constant fudge factor in an effective three- 
level atom. We find that the best fitting fudge factor 
would be F = 1.126, with relative errors reaching 0.2 
%. In any case, the effective MLA computation is so 
simple and computationally efficient that the need for 
non-physical fudge factors does not arise. 

In Fig. [4] we show the effect of adding higher order 
Lyman transitions and feedback between them. It is suf- 
ficient to include Lyman transitions up to Ly7, as ne- 
glecting higher transitions leads to relative changes of 
10~ 5 only [15 . This initially speeds up recombination 
by adding more decay paths to the ground state, then 
slows it down due to delayed reabsorptions of Ly/3 pho- 
tons in the Lya line. Our results are similar to those of 
Ref. [IT]. 

For now on, our "base" model will be the effective 
multi-level atom model with Lyman-a, /3 and 7 transi- 
tions and feedback between them (assuming a blackbody 
radiation field incoming on Ly7). In the next sections we 
graft two-photon processes onto this base model. 



IV. TWO-PHOTON PROCESSES: FORMAL 
DESCRIPTION 

A. Overview 

It is well known since the first works on primordial 
recombination that 2s — > Is two-photon decays signifi- 
cantly contribute to the recombination dynamics [TJ [2] . 
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FIG. 3. Fractional changes in the ionization history relative 
to the effective three-level atom model. The "RecFast" model 
is an effective three-level atom with the case-B recombination 
coefficient multiplied by a fudge factor F = 1.14. The same 
prescription for the evolution of the matter temperature is 
used in all cases, see SecIV El 
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FIG. 4. Fractional changes in the ionization history when in- 
cluding higher-order Lyman transitions and feedback between 
them, compared to the effective multi-level atom model with 
2 s and 2p only. 



Even with a relatively low decay rate, the forbidden 
2s — » Is decays are indeed comparable in efficiency to the 
highly self-absorbed Lyman-a transition; in fact, more 
than half of hydrogen atoms have formed through the 
2s — » Is channel jJD]. This process was traditionally ac- 
counted for with the total 2s — > Is decay rate in vacuum, 
A-2 S ,is ~ 8.22 s _1 , with a two-photon absorption rate ob- 
tained by detailed balance considerations. For the level 
of accuracy required for future CMB experiments, one 
needs to account for stimulated two-photon decays [22] 
and non-thermal absorptions 23, 24 . 

Recently, it was suggested that two-photon decays 
from higher lying ns and nd states may also lead to per- 
cent level corrections to the recombination history |25j . 
Inclusion of such decays presents an additional concep- 
tual difficulty which was not present for the 2s — > Is 
decays: the problem of double-counting. Indeed, there 
is no fundamental difference between a sequence of two 
allowed one-photon transitions nl — » n'p, n'p — > Is, with 
1 < n' < n, and a two photon decay from the nl state 
near resonance (i.e. where the energy of the two photons 
are near E nn > and E n i\ respectively). Approximate so- 
lutions were presented in Refs. [55J HI] (for a review, 



see Ref. [53]). The double-counting problem as well as 
the reabsorption problem were resolved with a numeri- 
cal approach, solving the radiative transfer equations for 
the photon field in Ref. [UJ , which also provided analytic 
approximations to check the validity of the numerical re- 
sult. In this work, we will use the same numerical method 
as in Ref. 24j, which we then extend to account for fre- 
quency diffusion near the Lya line. In this section, we 
review the formalism presented in Ref. [21] and how to 
solve the double counting problem. In Sec. [V] we will de- 
scribe our numerical method for solving simultaneously 
the radiative transfer equation and the evolution of the 
atomic level populations. 



B. Two-photon decays and Raman scattering 



We start by defining the coefficient: 



dv 



cef s v 3 v' 3 
108(2/ + Ipjf 



\M(v)f 



(47) 



where the matrix element M.{v) is given by Eq. (B5) of 
Ref. |24j . and v' = \v — v n \\, where v n \ is the frequency 
of the Ly-n transition. For v < i> n i, dA n i/dv is the rate 
of spontaneous two-photon decays from nl per frequency 
interval. For v > v n \, dA„i/dv x f v i (where /„ is the 
photon occupation number at frequency v) is the rate 
of spontaneous Raman scatterings per frequency inter- 
val per atom initially in nl (in the notation of Ref. [24] , 
dA n i/di> — dK n i/dv for v > v n \). The function dh n i/dv 
is continuous across v = v n \, where it vanishes. 

We can now write the net rate of nl O Is two-photon 
transitions per frequency interval per hydrogen atom, for 
which the highest energy photon has frequency v < v n \. 



dA 7 



dv 

gis 



Znl(l + /„')(! + /") 



(48) 



For v > i> n \, the appropriate rate is that of Raman scat- 
tering events: 



Ani(^ > V n \) 



dh-nl 

dv 



9u 



Xnlfu>(l + fu) 
-Xi s (l + f v >)fv 



(49) 



In both cases, we can assume that the photon occupation 
number for the low-energy photons is that of a blackbody, 
since the optical depth for two-photon absorption of the 
low-energy photons is tiny (for a discussion, see Ref. [24]). 
We therefore set /„/ = (e hv / Tr — l) -1 . Moreover, the 
photon occupation number for frequencies v > v^ ya /2 is 
much smaller than unity: /„ « 1. This means that we 
can neglect stimulated emission by the high-energy pho- 
tons in Eqs. (48) and (49 1. Given these considerations, 



the net rate of two-photon transitions can be written in 
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the following form, valid for both v < v n \ and v > v n \\ 



Near Lyman resonances, we have 



A ni (i/) 



dL 



dv 



nl\ e h(v-u nl )/T, 



r , _ M Tl P fe(f-^i)/r r 
31s 



(50) 



dv 



(v w i/ nl ) w 3A„ Pils p™ ^„(i/), 



(55) 



C. Resonant scattering in Lyman-a 



We now consider pure scattering events, 



H(ls)+ 7 ^H(ls)+ 7 . 



(51) 



In the low-frequency limit this corresponds to the famil- 
iar Rayleigh scattering phenomenon; the cross section 
however has resonances at the Lyman-series lines, which 
correspond to resonant Rayleigh scattering. 

Rayleigh scattering events conserve the photon fre- 
quency in the atom's rest frame. In the comoving frame 
(frame in which the CMB appears isotropic), the fre- 
quency of the scattered photon appears shifted due to 
the thermal motions of the scatterers. The frequencies 
of the incoming and outgoing photons are however sta- 
tistically correlated. Mathematically, there is a definite 
probability distribution p{v,v'), such that p(v,v')dv' is 
the probability that the outgoing photon has frequency 
in [z/, v' + dv'] given that the incoming photon had fre- 
quency v, and this function generally depends on both v 
and v' . For T m <C hv, which is the case near the Lyman 
lines, the variance of the frequency shifts imparted by 
thermally moving atoms is given by: 



(5v 2 ) = {v' - v) 2 p{v, v')dv' 



2T, 



m v 2 . 



(52) 



The rate of injection of photons per frequency interval at 
frequency v, due to resonant Rayleigh scattering in Ly-a, 
can be written in the general form (neglecting stimulated 
scatterings) : 

A la (i/) =x ls [f f v ,R{v', v)dv' - J f v R{v, v')dv'} , 

(53) 

where R(v, v') — ^j7f-p{v, v') is the differential rate of 
scatterings per hydrogen atom in the ground state, per 
unit frequency interval for both the incoming and out- 
going photons (it has units of s _1 Hz~ 2 ). The scattering 
kernel must respect detailed balance: 



R{v, z/)e-'W T '" = R{y' , v)e~ hv ' ' T ™ 



(54) 



To be fully general one should compute the scattering 
kernel from first principles. However, simplifications can 
be easily made in various regimes. 

Far from any resonance, the rate of redshifting due 
to the Hubble expansion is much larger than the rate 
of frequency diffusion due to scattering (see for example 
the discussion in Ref. [H]). We can neglect Rayleigh 
scattering there, and set Ai s (v) = 0. 



where p" c = A np ls /T np is the scattering probability in 
the Lyman-n line (the complementary events being two- 
photon absorptions and two-photon photoionizations), 
and 4>v,n{v) is the Voigt profile for the Ly-n line. In the 
Doppler core, we can approximate the partial redistribu- 
tion induced by scattering events by a complete redistri- 
bution, i.e. approximate p(v,v') §v,n(y) ~ 4>D,n{v), 
where (f>jj is the Doppler profile. This approximation 
is valid because in the Doppler core, complete redistri- 
bution recovers the correct rms frequency shift during 
scattering events, Eq. ( [52] ) (if one averages over the fre- 
quencies of absorbed photons). 

In the damping wings of Lyman resonances above Lya, 
the rate of scatterings is of the same order as the rate of 
two-photon absorptions. Each scattering event shifts the 
photon frequency by a very small amount compared to 
the width over which the radiation field varies (5v tms /v ~ 
2.5 x 10~ 5 ). Partial redistribution is therefore essentially 
coherent in the comoving frame, i.e. p{v, v') rj 8{v' — 
v), which implies Ai s (u) ~ 0. For a more quantitative 
argument, see Ref. [T5"] . 

The only frequency regime where Rayleigh scattering 
affects the radiation field in a non-trivial way is in the 
damping wings of Lya. In this line, indeed, scattering 
events are much more frequent than two-photon absorp- 
tion events (by a factor of ~ 10 4 ). Resonant scatter- 
ing therefore leads to a significant diffusion in frequency. 
Because the frequency shifts are small compared to the 
width over which the radiation field varies, the integral 
scattering operator can be approximated by a second 
order differential operator - a Fokker-Planck operator 
[2171 |2"T1 l4"2l H3] . For the purpose of numerical implemen- 
tation, the relevant properties are (i) the fact that this 
operator is nearly local (it only connects neighboring bins 
in frequency) (ii) it must respect detailed balance and 
(iii) the diffusion rate must be correct. We will explain 
our numerical method for the implentation of Lyman-a 
diffusion in Sec. IV Al 

We note that a number of analytic treatments of 
Lyman-a scattering in the recombination epoch have 
been proposed in the past [41&I47| . However, since two- 
photon emission and absorption act on the same region 
of frequency space, and since both processes involve high 
optical depth, an accurate recombination history can 
only be obtained by considering all processes simulta- 
neously. 
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D. The radiative transfer equation 

The radiative transfer equation for the photon occupa- 
tion number is: 



dfu 

dt 



Hv 



df„ = 

dv 8irv 2 



c 3 n H 



&nl{v) + Ai fl (l/) 



n>2,l 



, (56) 



where the left-hand-side is the derivative of the photon 
occupation number along a photon trajectory in the ex- 
panding universe, and the prefactor on the right-hand- 
side converts the number of photons per unit frequency 
per hydrogen atom to the photon occupation number. 



E. Inclusion in the effective multi-level atom rate 
equations 

1. Formal two-photon decay rates 

As mentioned earlier, including two-photon decays 
from states with n > 2 and Raman scattering events 
poses a double-counting problem. In principle, to avoid 
this double counting issue, one should discard "1+1" de- 
cays (or decays following an absorption event, which is 
like a Raman scattering event on resonance) altogether. 
If one were to pursue this idea, one should not consider 
the p states at all anymore (as they are formally only 
intermediate states in two-photon processes), but con- 
sider all s and d states as "interface states" and allow 
for two-photon recombinations to the ground state. The 
two-photon nl o Is transition rates would then become: 



I (27) 



-Xls 



I (27) 



— t, Rtotal „ , ntota 



total 



(57) 



where the formal transition rates are given by: 



ototal 



_ f dA n i g nl , h(u nl -v)/l 

dv g ls 



l| f v dv (58) 



and 



R 



total 
nlAs 



dk 



nl\ e h(u-v nl )/T T 



1 



dv, 



(59) 



where the integrals run from v n \/2 to v c . In principle 
Eq. ( 57 )— ( 59 1 , can be included in a standard or effec- 



tive multi-level atom code, provided one solves simulta- 
neously for the radiation field, using the radiative transfer 
equation Eq. (56). 



Decomposition into "1 + 1" transitions and non-resonant 
contributions 



Two-photon decays from higher excited states consti- 
tute, however, a correction to the recombination history 
computed in the standard "1+1" picture, and we would 



like to implement it as such. We start by formally sep- 



arating the integrals in Eqs. (58) and (59) in two con- 
tributions: the resonant pieces, for v s» iVii an d a non- 
resonant piece, for frequencies far enough from any reso- 
nance. We therefore rewrite, formally: 



R\°s,nl = X! ^U,nl + Ru,nl and 



ptotal _ f>(n'p) , p 
n nl,ls — 2-^i n nl,U + n nl,U 



(60) 



n' 



where the resonant contributions ^) and R^ are 



defined in a similar manner as in Eqs. (58) and (59), but 



with the integration being carried over a narrow range 
Av near v n 'i, and Ri s , n i and R n i.\s are the non-resonant 
pieces required to complete the total rates. So far the 
separation is just formal and we have not made any ap- 
proximation. 



3. "1+1" Resonant contribution 

We now notice that near a resonance v w ivi, the 
two-photon differential decay rate dA n i/dv takes on the 
following form (if n > n'): 



dA„ 



dv 



A 



nl.n'p-^n'p^ls 



An 2 (v - ivO 2 + (W4tt) 2 



P,ls 



4>l{v - v n n;T n > p ), (61) 



where T n i p is the total inverse lifetime of the state n'p, 
and the Lorentzian profile is given by 



i z (Ai/;r) = 



r/(4^ 2 ) 

Av 2 + (r/47r 



|2 ' 



(62) 



For n < n', the first coefficient in Eq. (61) should be 



gn' P /gni x A n ,p^i instead of A n i t n> p . When accounting 
for the thermal motions of atoms, the Lorentzian pro- 
file should be replaced by a Voigt profile. We can now 
approximate the resonant pieces with the following ex- 
pressions, valid for both n < n' and n > n': 



r>( n 'p) OA 
"■13,711 ~ OJ± n'p,lsJ Unll 



R 



n'p,nl 



r, 



and 



R 



(n'p) 
nlAs 



RnLn' 



x n'p 



(63) 



(64) 



where /„ is the photon occupation number averaged 
over the Voigt profile near the resonance v w v n '\- 
Eqs. (63) and (64) are exactly what one would obtain in 



the "1+1" picture after "factoring out" the p states (with 
a procedure similar to what is used to get rid of the "in- 
terior" states in the EMLA method). Having these res- 
onant rates is exactly equivalent to having optically thin 
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one-photon transitions between the nl and n'p states, 
with rates R n i,n'p(T T ) and Rn'p,ni(T T ), and optically thick 
Lyman transitions, with net rate: 



contributions of two-photon decays is accurate to better 
than 10~ 4 . 



bn p\u 



\%xuf Vn , p - x n ,j)j . (65) 



To obtain the net decay rates in the Lyman transitions, 
one then needs to solve for the radiation field in the im- 
mediate vicinity of Lyman resonances. If the frequency 
region for which two-photon transitions are considered 
as "resonant" is narrow enough, this can be done in the 
Sobolev approximation. Indeed, all the relevant condi- 
tions are met (see also discussion in Ref. [23]; for more 
details on the Sobolev approximation, see for example 
Ref. [H]): 

First, the two-photon absorption and emission profiles 
can both be approximated by the same resonance profile 
Eq. (61). This relies on the assumption that the black- 



body radiation field varies little across the "resonant" 
region, and requires for its width to satisfy Av <C T T /h. 

Secondly, we argued in Sec. |IV C] that one could assume 
complete frequency redistribution for resonant scattering 
near the Doppler core of Lyman resonances. The "reso- 
nant" region should therefore not exceed a few Doppler 
widths. 

Finally, if we consider regions in frequency narrow 
enough around the resonances, we can use the steady- 
state approximation. This requires Avjv <C 1. 

We can see that considering the "resonant" region 
around each Lyman resonance to be a few Doppler widths 
wide meets all the requirements. 

An additional assumption required here is that excited 
states are near Boltzmann equilibrium, which is very ac- 
curate at redshifts for which two-photon processes are im- 
portant. In the Sobolev approximation, and in the limit 
of large Sobolev optical depth, Eq. (|65l) becomes the stan- 



dard Lyman decay rate Eq. (|18|), where /+, is the photon 



occupation number incoming on the resonance, prepro- 
cessed by two-photon processes and diffusion in the blue 
damping wing of the line. 

The Sobolev approximation is probably the least ac- 
curate for Lya decays where partial redistribution due 
to resonant scattering is important. However, the large 
optical depth to two-photon absorptions in the Lyman-a 
blue damping wing, in conjunction with frequent scat- 
terings, drive the radiation field to the equilibrium value 
fu = x n i p /{ixi s )e~ h ^~ Un ' 1 ^ TTn over several Doppler 
widths (of the order of 40 Doppler widths, see Ref. [15]). 
As a consequence the net decay rate in the core of the 
resonance is very small anyway. We checked that in the 
presence of two-photon transitions and frequency diffu- 
sion, even setting i„' p | ls = instead of the expression 



given by Eq. ( 65 ) leads to relative changes to the recombi- 
nation history of at most 7 x 10 -4 . Given that frequency 
diffsion leads to corrections of a few percent at most to 
the decay rate in Lya when radiative transfer is treated 
carefully even at the line center [2U] , we can be confident 
that using the Sobolev approximation for the resonant 



4- "Pure two-photon" non-resonant contribution 



In the previous section we discussed how two-photon 
decays within a few Doppler widths of Lyman resonances 
can in fact be accounted for in the standard "1+1" pic- 
ture. To evaluate the non-resonant pieces, Ri s , n i and 
Rni,is, we need to solve the radiative transfer equation, 



Eq. (561, to obtain the photon occupation number. The 



subject of Sec. [V] is to describe our numerical method of 
solution. 

Note that choosing the "resonant" regions to be a few 
Doppler widths has an additional advantage. Since a 
Doppler width is ~ 10 3 times wider than the natural 
width of Lyman lines, it is not necessary to account for 
the pole displacements in the computation of the differ- 
ential two-photon decay rates in the non-resonant region. 
In addition, the fraction of two-photon decays that are 
considered non-resonant will be small (of the order of 
r np /(47r 2 )/Az/, where Av is the width of the "resonant" 
region). For Av of a few Doppler widths, this fraction 
is ~ 1CP 4 . This means that the "pure" two-photon de- 
cay rates R n i,i s are much smaller than the total inverse 
lifetime of the nl state, r„;, which is required to simplify 
the effective MLA model to an effective four-level atom 
model as we discussed in Sec lIIIB"2l 

As a final note, we want to emphasize why the final re- 
sult is independent of the exact boundary between "res- 
onant" and "non-resonant" regions, so long as the res- 
onant regions are a few Doppler widths wide. If one 
were to increase the width of the "resonant" region, then 
the "pure" two-photon transition rates R n i,i s and Ri s<n i 
would decrease, mainly because of the change of the inte- 
gration region in the blue wings of the resonance - in the 
red wing, the radiation field has reached near equilibrium 
with the line and the net rate of decays immediately blue- 
ward of line center is very small anyway. This decrease 
would be nearly exactly compensated by the increase of 
what is considered as "1+1" decays, as the photon occu- 
pation number incoming on the Lyman resonances, 
would be decreased due to the smaller optical depth due 
to "pure" two-photon absorptions in the blue wing. Hi- 
rata (2008) checked the independence of the result form 
the exact value chosen for the width of the "resonant" 
region, and found that even changing this width by a 
factor of 9 lead to relative changes of at most 4 x 10~ 4 
in the recombination history. 
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V. NUMERICAL SOLUTION OF THE 
RADIATIVE TRANSFER EQUATION 



these conditions at the boundaries of the diffusion do- 
main, smaller than the entire frequency domain consid- 
ered). Using Eq. (70), we then obtain Ri 2 an d Rn,n-i- 



A. Discretization of the radiative transfer equation Using iteratively Eqs. (69) and (70), we can then obtain 



To solve the radiative transfer equation [Eq. (56)] nu- 
merically in the "non-resonant" frequency region, we fol- 
low the method of Hirata (2008), and extend it to also 
account for frequency diffusion. 

We will consider the radiation field in the vicinity of 
N frequency "spikes" v b , for b = 1,2,... N. Each spike 
has an associated width Avb (which is just the separation 
between consecutive spikes if they are linearly spaced for 
example) . 

We use the discretized differential two-photon rate: 



dA n i 



dv 



= 22 A rdjb 5 e (v - v b ), 

used — 
b 



(66) 



where we use the coefficients 



all the coefficients of the numerical diffusion kernel, start- 
ing from the boundaries, and up to line center. Denoting 
b\ the highest bin below Ly-a and 6^ + 1 the first bin 
above Ly-a, we obtain all coefficients up to Rb^ ya ,b x on 
the red side of Ly-a, and up to -R& LyQ .6i+i on the blue 
side (we do not follow the radiation field at the central 
bin &Lya but can still define these coefficients). Note that 
with this method we cannot ensure that the diffusion rate 
at the central bin is correct. However, the exact value of 
the diffusion rate at line center does not matter, as long 
as it is high enough to ensure that the photon occupation 
number reaches the equilibrium spectrum f v oc e~ tlv / T <^ . 



B. Solution of the discretized radiative transfer 
equation 



-4, 



Lb 



dA r , 



-dv, 



(67) 



where the integral is carried over the frequency region 
associated with the spike Avb- The function 5 t {v — z/&) in 
Eq. ( 66 1 should be understood as a sharp profile centered 
at Vb, which integrates to unity, and has support in [vb — 
e, Vb + e]. The solution we derive is in the limit e — > 0, for 
which 6 e — > S, the Dirac delta function. Doing so, we are 
simply approximating the optical depth as concentrated 
in discrete frequencies instead of being a smooth function. 

The main new contribution of the present work is the 
discretization method for the scattering operator. We 
use the discretized scattering kernel 



R(v, v'] 



used 



E 

6,6' 



Rb,b' $e {v - v b )S e {v - v v ) ■ (68) 



We enforce detailed balance: 



Rb,b'C 



-hv b /T m _ 



Rb'.bZ 



(69) 



We moreover use the diffusion approximation for reso- 
nant scattering. This allows us to assume that the nu- 
merical scattering kernel Rb.b' is non-vanishing only for 
neighboring bins, b' = b±l. In order to obtain the correct 
diffusion rate, we set 



(yb+i 
= 3 



Vb) 2 Rb,b+l + iyb-\ — Vb) R b ,b-1 



A 2 



2T„ 



Aw 2 (v - v hya y 



A— in z 
Vb ^Lya, 



TOfjC' 



(70) 



where we used the damping wing approximation for the 
absorption profile (and approximate v 2 m v^ ya in the 
multiplicative factor). 

As boundary conditions, we assume a vanishing pho- 
ton flux due to diffusion at the boundaries of our do- 
main, i.e., formally, Ri t o = Rn.n+i = (in fact we set 



To simplify the notation, we define the following rate 
coefficients: 



R, 



dA r , 



Lb 



Ji{v b -v nl ) /T r 



Rb,: 



dv 



9nl c h{ Vb - Vnl )/T r 

gis 



1 Avb and 



Ri 



Lb- 



(71) 



This coefficients can be thought of as transition rates 
between bound states and a set of "virtual" levels with 
associated energies Eb = hvb [2~i] . 

We define the total Sobolev optical depth in the &-th 
frequency spike: 



SirvgH 



(72) 



nl 



b'=b±l 



We also define the average photon occupation number 
near v^. 



f„ 



!'):, + <■ 



5 e (v - v h )f v dv. 



(73) 



Finally, we define the equilibrium photon occupation 
number at the 6-th frequency spike: 



feq 

J Vh 



J2nl x nlRnl,b + Xls J2b' f v, , ^6' 



6,6' 



(74) 



In the vicinity of Vb, the discretized radiative transfer 
equation becomes: 



1 df v df v 



Hvb dt dh 



Ar b 5 e (v-v b ) [Q-U] 



(75) 



In the limit that the support of the delta function be- 
comes vanishingly small, e — > 0, the discretized radiative 
transfer equation can be solved in the steady-state ap- 
proximation, and one can neglect the time derivative. 
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This is similar to the commonly used Sobolev approxi- 
mation, except that we are now making this approxima- 
tion in the vicinity of an artificially introduced spike (as 
opposed to a true resonance line), for the purposes of nu- 
merical resolution. Another conceptual difference is that 
the equilibrium photon occupation number also depends 
on the averaged value of the radiation field at neighboring 
bins, because of frequency diffusion. Given the photon 
occupation number at the blue edge of the 6-th spike, 
fv b +e> this equation has a well known solution f v . The 
quantities of interest for us are the photon occupation 
number at the red edge of the spike /„ b _ e and the aver- 
age photon occupation number in the spike f v . They 
are given by the following expressions (for a derivation, 
see for example Refs. [T5] and |24) ) : 



fvb — f- /' 



and 



f Vb = n b u b+e + (i - lit)/" 



(76) 



(77) 



where Ilf, is the Sobolev escape probability from the fr-th 
spike: 



IT 



1 - c- ATb 
An ' 



We now use the variables 



Xb = x u f v 



(78) 



(79) 



As explained in Ref. [21], Xb can be interpreted as the 
population of the virtual level b. One should however 
keep in mind that this is is just a convenient rewording 
for the radiation field intensity. 



Using the definition of /° q , Eq. (74), we can rewrite 



Eq. ( 77 ) in the form: 



T b ,bXb = ^2 X nl R nl,b + ^2 XyRb'.b + S b , (80) 



b'=b±l 



where we have defined: 



Y~n h ( 12 Rb > nl + 12 Rbb ' ) 

° \nl b'=b±l / 



Tb,b 

Sb = TlbXisfv b +cT b ,b 



and 



(81) 



We only follow two-photon decays in the damping 
wings of resonances, but we should still account for fre- 
quency diffusion between line center and the neighboring 
bins. At the Lyman-a line center, the radiation field is 
in equilibrium with the 2p/ls ratio: = X2 P /(3xi s ). 

If b\ is the highest frequency bin below Ly-a, (and b\ + 1 
is the first bin above Ly-a), we therefore define the tran- 
sition rates with 2p: 



R 



1 



2p,6i 



~Rb 



1 



<bl and J?2p,6i+1 = q^ fc Ly Q ,6i+i- (82) 



Provided that we set Rb 1 ,b 1 +i 



Rb 



-l,6i 



0, Eq. (80) 



remains valid for b = b\, b\ + 1. Adding these transitions 



with the central frequency bin will ensure that the photon 
occupation number is driven to its equilibrium value near 
line center, = W(^i s )e-' l(l '"'' Ly ° )/r °. 

We now use Eq. ( [42] ) for the populations x n i with n > 
3. We define the coefficients, for i = 2s, 2p: 



T b , = -Ri, b - J2 —e- E - 2/Tr P^(T r )R nlib . (83) 



9nl r -E n2 /T t pi 

n>3,l 

We define the new source vector: 

S b = s b + x\ J2 9 —e- E -^P e nl {T,)R nl ^ (84) 

n>3,i ^ e 

We also define the coefficients 

T b ,b±l = —Rb±i,b- (85) 

The discretized radiative transfer equation then takes the 
final form: 



6+1 

E 

b'=b-l 



T b ,2sX2s + T bt 2 p x 2p + y2 T b,b'X b ' = S b - (86) 



C. Populations of the excited states 



Given the radiation field, we can now compute the two- 



photon transition rates. Using Eqs. (58) and (59) limited 



to the "non-resonant" frequency region, we o 
discretization: 



RnlA 



12 R n 



Dtain, after 

/./, and Ru,nl = f Vh R b,nl- (87) 



The effective transition rates from the i = 2s and 2p 
states to the ground state are therefore, according to the 



discussion in Sec. Ill B 2 and using the definition of T b ^ 
Eq. ([831: 



Ki,U = -12 T b ,i + J2 9j ^ R Lyne- E ^PU T r), 



where the first term accounts for two-photon transitions 
and the second term for escape from Lyman lines (it is 
understood that = and P^ = 1). The effective 
transition rate for the reverse process is given by: 

n ls ,i = -J2 Ti,bf„ b + 3 ]T R Lyn Pi p (T I )f+, (89) 



where we have defined the coefficients, for i = 2s, 2p: 



T, 



Lb 



9i e h(v b —u 2 -i)/T r rp 



9la 



b.i- 



(90) 



For a given radiation field and free electron fraction, we 
can now obtain an equation for the populations of the 
excted states 2s, 2p. We do so in the steady-state ap- 
proximation, i.e. setting ii — for i = 2s, 2p in Eqs. (37) 



and ( 38 ) . We first define the 2x2 matrix of elements 
T;i = B t + TZii + Hiis and 



(91) 
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FIG. 5. Sparsity pattern of the linear system solved for evolv- 
ing simultaneously the level populations and the radiation 
field, in the presence of two-photon transtions and frequency 
diffusion. 



We also define the source vector of elements 

Si ee n n xlA + 3 £ iW^(T r )/+ . (92) 



n>2 



The steady-state equation for each state i translates into 
the linear equation: 

^ TijXj + T ifi x b = Si, i = 2s, 2p. (93) 

j=2s,2p b 



D. Evolution of the coupled system of level 
populations and radiation field 

We now have all the necessary pieces to evolve simul- 
taneously the level populations and the radiation field, 
and compute the free electron fraction. In this section 
we summarize the procedure and recall the main equa- 
tions. We start with an initially thermal radiation field, 
f u = or hv l Tl . At each time step, we do the following 
computations: 

1. We obtain the photon occupation number incom- 
ing on each bin b assuming free streaming between 
frequency spikes: 



fv b +e(z) — f, 



fb+i- 



(l + z) 



Ei 



6+1 



E b 



1 



(94) 



We also obtain in the same way the incoming pho- 
ton occupation number at the Ly-n transitions, 

2. We solve for the populations of the 2s and 2p states 
and the average photon occupation number at each 
frequency spike /„ = Xb/x\ s simultaneously by 
solving the coupled linear system given by Eqs. ( 86 ) 
and ( 93 ) . Even with a large number of bins for the 
radiation field (N = 311 in our fiducial case), this 



system is easily solved because the matrix of coef- 
ficients Tf)^ 1 is triadiagonal and the overall system 
has the particular sparsity pattern shown in Fig. [5] 
Such a sparse system can be solved in O(N) op- 
erations (specifically, we can solve the system in 
~ 16iV operations). 

We update the photon occupation number at the 
red side of each spike, f Vb - e , using Eq. (76). At 



the red side of Lyman resonances, we use /, 



np 



j/(3xi s ), valid in the optically thick limit, where 



is given by Eq. ( 42 ) for n > 3 



After step =#p] we can ob tain the function x e (z, x e ) 
through Eq. ( f39| oi[^] ( |40| . This allows us to evolve 
the free electron fraction to the next timestep. 



E. Implementation, convergence tests and results 

We evolve the free electron fraction during hydrogen 
recombination in several phases. We use even steps in In a 
(where a = (1 + z) _1 is the scale factor), with Alna = 
8.5 x 10 -5 . We describe our ODE integrator in Appendix 



C. 



• We checked that hydrogen and helium recombina- 
tion never overlap and can be followed separately 
(to an accuracy of a few times 10 -4 ). We therefore 
start computing the hydrogen recombination his- 
tory once helium is completely recombined. Quan- 
titatively, we start hydrogen recombination once 
the fractional abundance of He + ions is less than 
10~ 4 relative to hydrogen. If this criterion is met 
earlier than z = 1650 we only switch on the hydro- 
gen recombination computation at z — 1650. We 
checked that at this redshift the exact free electron 
fraction differs from the Saha equilibrium value by 
no more than a few times 10~ 4 anyway. 

• In the first phases of hydrogen recombination, we 
use the post-Saha expansion described in Appendix 
|D| We do so as long as the free electron fraction 
differs from the Saha value by Ax e < 5 x 10~ 5 . We 
checked that explicitly integrating the ODE for x e 
instead (with a much smaller timestep as the ODE 
is stiff at early times) leads to maximum changes 
of Ax e /x e < 3 x 10" 4 . 

• From then on and until z — 700 we solve simul- 
taneously for the level populations and radiative 
transfer with two photon processes and diffusion as 
described in this section. 



2 Eq. j39| l contains near exact cancellations but Eq. (|40| contains 
a large number of terms for which numerical roundoff errors can 
add up. We checked that both equations give the same result 
within numerical roundoff errors. We use Eq. ( |39| in the final 
code simply because it is more compact. 
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• For z < 700, we use the simple EMLA equations, 
with simple decay rates from 2s and 2p only (i.e. 
not accounting for higher order Lyman lines and 
radiative transfer effects) . Wc checked that moving 
the last switch to z = 400 instead of 700 leads to 
maximum changes |Ax e |/x e < 10 -4 . 

• For the matter temperature evolution, we use the 
asymptotic solution of Hirata [21] [it can be ob- 
tained by setting T m = —HT m in Eq. (12)] as long 
as 1 —T m /T T < 5 x 10~ 4 . Depending on cosmology, 
this corresponds to 750 < z < 950. After that we 
switch to solving for x e and T m simultaneously by 
using Eq. (fl2l. 



All the checks mentioned above were made for a wide 
range of cosmological parameters. 

Our fiducial parameters for the numerical solution of 
radiative transfer are N = 311 frequency bins extending 
from ^Lya/2 to Vhyy, and a diffusion region with 80 bins 
extending to Av jv^ya = ±1.7 x 10~ 2 . The minimal spac- 
ing between bins is min[ln(i/ft + i/ft)] = 8.5 x 10~ 5 , which 
sets the largest step in In a that we can take (this is also 
half of the width Av/v of the "resonant" region around 
Lya). We checked (for the fiducial cosmology only) that 
reducing the diffusion region to Avjvx, ra = ±1 x 10 -2 
leads to changes \Ax e \/x e < 6 x 10 -6 . Reducing the dif- 
fusion region to Ai//^Lya = ±5 x 10 -3 leads to changes 
\Ax e \/x e < 4 x 10~ 5 . We checked that using a 10 
times finer frequency grid in the diffusion region (and 
a 10 times smaller timestep) leads to maximum changes 
\Ax e \/x e « 1.5 x 10~ 4 at z « 900. 

We are therefore confident that our numerical treat- 
ment is converged at the level of a few parts in 10 4 . 

We show in Fig. [6] the changes in the recombination 
history due to two-photon processes. We find that in- 
cluding two-photon transitions from the initial states 
2s, 3s, 3d, 4s and 4d is sufficient for the level of accu- 
racy required - we checked that including two-photon 
transitions from 5s and 5d leads to a maximum change 



Ax e / x e 



8 x 10- b at z 



1200 and can therefore be 



neglected. The effect of frequency diffusion in the Lya 
line is shown in Fig. [7j 

We compared our results to the two-photon MLA code 
of Hirata [24] , as well as to the results of Hirata & Forbes 
[2"0] for frequency diffusion. For this comparison, we use 
n max = 30. The result of the comparison is shown in 
Fig. [8] The maximum difference between the codes for 
700 < z < 1600 is \Ax e \/x e = 0.0005. The increase 
of the relative difference at late times is most likely due 
to small differences in the bound-free rates, which are 
computed with different methods (we use the recursion 
relations of Ref. [48] whereas Hirata (2008) directly inte- 
grates the products of wave functions to compute matrix 
elements). This difference remains even when switch- 
ing off two-photon processes. The ~ 3 x 10~ 5 kink at 
z ~ 1570 is a startup transient due to switching from the 
post-Saha solution to solving the full ODE. The kink at 
z = 1350 is due to the Ly/3 photons emitted at z = 1600 
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FIG. 6. Changes in the recombination history when includ- 
ing two-photon decays and Raman scattering (no diffusion), 
compared to our "base" model. The line labeled '2s' shows 
the changes in x e when one properly accounts for stimulated 
2s — > Is decays, as well as absorptions of distortion pho- 
tons and Raman scattering form 2s. The other lines show 
the cumlative correction when adding two-photon transitions 
from higher levels. 




1600 



FIG. 7. Changes in the recombination history when including 
frequency diffusion in Lyman-a, compared to a model with 
two-photon transitions but no diffusion. 



starting to redshift into Lya. Overall, the agreement is 
excellent (\Ax e \/x e < 10~ 4 for z > 900), even though 
the codes use different methods to compute atomic rates 
and different approaches for solving the MLA problem 
and treating Lya frequency diffusion. 



VI. HELIUM RECOMBINATION 

Roughly 14% of the electrons in the Universe are asso- 
ciated with helium rather than hydrogen, so it is critical 
to model helium recombination as well [33 35 . Be- 
cause the ionization energy of helium is greater than that 
of hydrogen (E Il = 24.6 eV for He I and E l2 = 54.4 eV 
for He II) , helium recombines earlier than hydrogen, well 
before the epoch of last scattering. Therefore, we do not 
observe helium recombination directly: rather it affects 
the diffusion of photons at early times and hence controls 
the Silk damping length [5 ■ A faster helium recombina- 
tion (smaller x e ) leads to a longer photon mean free path 
and hence a larger damping length. The net effect is then 
to reduce the high-€ multipoles of the CMB temperature 
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ionization energy of helium and g e is given by Eq. Q but 
using the appropriate reduced mass of the electron-He III 
system. We first solve Eq. (95) for q: 



2s/ Hc 



1 + /He 




4s/ : 



He 



(1 + /He + S) 



(96) 



We then obtain the free electron fraction from x e = 1 
/He + <?■ 

This equation is used to obtain x e until q = XHeiii 
10~ 9 (which corresponds to z ~ 4000). 



B. He II^I recombination: Near-equilibrium stage 



FIG. 8. Comparison of HyRec with the MLA+two-photon 
code of Hirata (2008) and to the results of Hirata & Forbes 
(2009) who account for Lya diffusion. See text for comments. 



and polarization power spectra |33j . However, since most 
of the Silk damping occurs at later times, we do not need 
extraordinary accuracy in following helium recombina- 
tion: the corrections identified by Refs. [3TH55] . which 
changed x e by up to 3% at z « 1800, amounted to a 
r~> la correcton for Planck and ~ 8<r for a hypothetical 
cosmic variance limited experiment to I = 3000. There- 
fore, a helium recombination code accurate to ~ 0.3% 
should reduce any residual errors to the point of being 
negligible for Planck. Here we describe our "fast" helium 
recombination code. 

Throughout the helium recombination process, we may 
take a single temperature T = T m = T T . We make 
several other crude approximations, as detailed below. 
While analytically motivated, their quantitative justi- 
fication rests on the comparison to the "full" version 
of our calculations |33j . For the equilibrium calcula- 
tions, the He I level energies and are obtained from 
the NIST database (itself based on the compilation of 
Ref. [H]). The bound-bound Einstein coefficients are 
obtained from Ref. [50]: A[2 l P° - l x 5] = 1.7989 x 10 9 
s" 1 and A[2 3 P° - \ X S\ = 177.58 s" 1 . The bound-free 
rates are not required here since the excited levels of He 
I remain in equilibrium with the continuum (i.e. these 
rates may simply be taken to be "fast"). 



A. He III— >II recombination 

The He III— ^11 recombination has previously been 
found to be in Saha equilibrium to high accuracy [311 134] 
and is followed using the Saha equation: 



g(l + /Hc + g) 
/He - q 



5 e e 



-EjJT 



(95) 



where q = XHciii and we assume that the rest of the 
helium is singly ionized, XHcii = /hc — 9j and all hydrogen 
is ionized, x e — 1 + /hc + <7- In Eq. (95), Ej 2 is the second 



In contrast to He III— >Il recombination, the He II ^1 
recombination is a highly non-equilibrium process (it oc- 
curs at much lower density and in a weaker radiation 
field, with a much slower 27 decay process, and the ex- 
cited levels of He I are much closer to the continuum as 
a fraction of the ionization energy than those of He II) . 
Therefore, it must be followed in several stages. 

The first is the near-equilibrium stage, before the main 
He I resonance line (2 1 P°-1 1 S at 584 A) develops suffi- 
cient optical depth to push helium recombination out of 
equilibrium. In this stage, the free electron fraction is 
close to the Saha solution x^, aha = 1 + q, where q = xudi 
satisfies 



/He 



J/) 

- 9 



(97) 



and we assume that all hydrogen is ionized and tha t th e 
helium is distributed between He and He + . In Eq. ( 97 1 , 
Ei 1 is the first ionization energy of Helium, and one 
should use the appropriate reduced mass of the electron- 
He II system in g e . The additional factor of 4 relative 
to Eq. ( 95 ) is due to the lower spin degeneracy of He I. 



Again, we first solve Eq. (97) for q 




(98) 



from which we get x^ = 1 + q. We then obtain the 
post-Saha expansion for the free electron fraction x e = 
x Saha _|_ ag described in Appendix [d] We do so until 
Ax e reaches 5 x 10 -4 , which corresponds to 2500 < z < 
3000 depending on cosmology. We checked that using the 
post-Saha expansion until Ax e = 10~ 5 instead and then 
numerically integrating the ODE for x e (described in the 
next section) leads to maximum changes of ~ 3 x 10~ 4 
in the free electron fraction. 



C. He II— >I recombination: Non-equilibrium stage 

At lower redshifts, one must follow helium recombina- 
tion carefully, via a Peebles-style ODE [I], i.e. an equa- 
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tion of the form 



1. The two-photon process 



J~(x e . 



(99) 



where p is the vector of cosmological parameters. Despite 
the complicated radiative transfer physics in the helium 
problem, a single ODE turns out to suffice because the 
portion of the ultraviolet spectrum that is relevant is rel- 
atively narrow and may be treated as in steady state (it 
encompasses the 584 and 591 A lines). The construction 
of the J- function is however more complicated than the 
Peebles ODE for hydrogen. 

The excited levels of He I (n > 2) have been found to 
remain in Saha equilibrium with the continuum through- 
out the process since the strong CMB blackbody ionizes 
them far faster than they can reach the ground state [3T] . 
The significant processes for net decays to the ground 
state are four-fold: 

• The two-photon process, He(2 1 S')-^He(l 1 5)+7+7. 

• Emission of photons via the main resonance line 
2 1 P°-l 1 S with A = 584 A, followed by redshifting 
out of the line (the He I analogue of the H I Lyman- 
a escape process). 

• The absorption of He I 584 A photons by H(ls) 
atoms via photoionization, H(ls)+7 — >H + + e~. 
The electron rapidly thermalizes its energy, leading 
to a loss of a resonance line photon and a net decay 
of He I. 

• Emission of photons in the intercombination line He 
I] 2 3 P°-1 1 S with A = 591 A. This line has Sobolev 
optical depth of order unity during He I recombina- 
tion so the full Sobolev escape probability formula 
must be used. 

In constructing the function J- above, we take several 
steps. First, we compute the abundance of the species 
H I, H II, He I and He II. The formulae for the latter 
are zhcII = x e - 1 + x m and a; Ho i = /tie - ^Hdi, but 
these require knowledge of the H I fraction. This cannot 
be completely igno red, but it is small and so the Saha 
equation for H (Eq. D2 ) suffices for this purposer] 

Next is the determination of the excited level (^S and 
2 1 P) abundances assuming equilibrium with the contin- 
uum, 



1 



-[215] 



i) 



(100) 



and x [2 i P] = 3x [2 i 5] exp [-(E [2 i P] - E [2 i s] )/T] . 

Armed with this information, we proceed to investigate 
each decay mechanism. 



3 Eq. IJD2JI technically neglects the contribution of helium to the 
electron abundance, but this correction is small in the very latest 
stages of He recombination when the H I correction becomes 
significant. 



The downward rate for the He I 27 decay is A =50.94 
s _1 [51] . so we find a downward decay rate 

(101) 



2. The 



line 



The rate of emission of photons in the 584 A line is 
complicated because photons in this line experience mul- 
tiple processes: (i) "true" absorption and emission, (ii) 
resonant scattering (which redistributes photons in fre- 
quency since in the comoving frame the ingoing and out- 
going photons need not have the same frequency), (iii) 
H I photoionization opacity (which has absorption and 
emission), and (iv) redshifting due to Hubble expansion. 
We construct here a highly simplified model of these pro- 
cesses; the "full" treatment can be found in Ref. |31| . 
Ref. [3T] also showed that the resonant scatterings can 
be neglected. It should be noted that in the case of the 
584 A line, the nontrivial physics takes place in the damp- 
ing wings: the line center is very optically thick, and the 
reaction 



He(2 1 F°) O He(l 1 S') +7 



(102) 



reaches equilibrium there (analogous to the H I Lyman-a 
line). 

In this situation, the radiative transfer equation for the 
photon occupation number f v in the vicinity of the 584 
A line can be simplified to [Eq. (36) of Ref. 31j, with 
individual expressions substituted in] 



dv 

-HvT ahs (j){v) (f^ 4 



(103) 



where rj c is the H I opacity in units of absorption optical 
depth per unit frequency (i.e. Hvr) c is the absorption op- 
tical depth per unit time), and f^ 4 = x\2i-p°]/{3x[iig]) 
is the equilibriu m pho ton occupation number in the 584 
A line. In Eq. (103 1, the e~ hv l T term corresponds to 
emission from direct H I recombinations to the ground 
state (i.e. the inverse process of photoionization); r a b s 
is the Sobolev optical depth to true absorption in the 
584 A line; 4>{v) is the 584 A line profile normalized by 
J 4>{v) dv = 1. We neglect stimulated emission processes 
at 584 A, since the photon phase space den sity is <C 1. 
In steady state, the left-hand side of Eq. (103) is ap- 



proximated as zero. The remaining contributions are as 
follows. 

The H I continuum optical depth is, assuming an H I 
abundance in Saha equilibrium with H II (assumed to be 
nearly all hydrogen, i.e. x p ~ 1), 



>k 



O-piCnuXe 



Hv 



(104) 
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[This comes from combining Eqs. (27) and (28) of leading to 
Ref. [31].] Here <r p i is the photoionization cross section 
of H I at A = 584 A and Ej is the hydrogen ionization 
energy. 

The line profile is in principle a Voigt profile. However where 
near line center where <j>{v) becomes large, we have f u — > 



^584 irrespective of the details in order to keep Eq. (|103[) 



finite. Therefore, we approximate it by the damping wing 
approximation, 



Att 2 (v - v Q ) 2 ' 



(105) 



where vq is the central frequency of the line and T is the 
intrinsic width. The latter is the sum of the rates for 
all processes that depopulate 2 1 P° (the width of the l 1 ^ 
state is negligible), i.e. we may write T = + r ot her, 
where A^4 is the contribution from 2 1 P° —> l 1 ^ and 
T ther is the contribution from all other states: 



.4 



other 



c -(E [2 i po] -E [2ls] )/T 

3 e (E,-E [2lpo] )/T _ 1 ' 



(106) 



where the first term corresponds to decays to 2 1 1S' (sup- 
plemented by stimulated transitions) and the second 
term corresponds to absorptions from 2 1 P° to a higher 
level i. We include the levels and n 1 !) for 3 < n < 5. 
In principle we should include higher n and the contin- 
uum levels, but in practice the first few levels dominate 
the sum because of the exponential factor. 

We next need the optical depth to true absorption, 
which is the Sobolev optical depth rs times the fraction 
of photon absorptions 1 X S — > 2 1 P° that are true ab- 
sorptions (i.e. do not immediately decay back to l 1 ^, 
but rather visit another level). This fraction is r ot h cr /r. 
Thus 



, / \ r other 
Tabs0(^) = ^ T S 



thcr 



(107) 

With these results, and neglecting the variat ion o f the 
blackbody function e ~ ,w / T across the line, Eq. ( 103 ) sim- 
plifies to 

9fu _ (f _ -hu /T\ , T S r othcr , , f cq, 

dv ~ Vc[tv 6 ) + 4^(v-v ) 2 [U hsi) - 

(108) 

The next simplification of this equation occurs if we re- 
scale both the frequency and the phase space density axes 



4tt 2 



r s r 



and 



othc 



-hvo/T 



f c.q 
Jl 



584 



-huo/T ' 



(109) 



(110) 



ay T 



TsTother^c 
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(111) 



(112) 



For blackbody radiation entering the line, we have the 
boundary condition £(+oo) = 0. 

This reduces the radiative transfer equation to a single 
ODE that depends on the single dimcnsionlcss parame- 
ter t c , which is (roughly speaking) the optical depth to 
H I photoionization within the part of the line that is 
optically thick to true absorption by He I 584 A. This 
parameter is exponentially increasing in the early parts 
of helium recombination, and becomes of order unity at 
z fa 2100. It is the parameter that controls which process 
is a more important sink for resonance line photons: H 
I continuum opacity (r c 3> 1) or escape via rcdshifting 

(r c « 1). 

The net rate at which photons are emitted in the He 
I line in photons per H nucleus per unit time is then 
obtained by integrating the absorption/emission term in 
the rate equation, 



jr(584) 



8tt^ 
n H c 3 



Hv T&hs( j>(v) (/ 5 c 8 q 4 - /„) dv. (113) 



This integral can be converted into an integral over y and 
f v can be replaced with £(y); making these substitutions 
and approximating v ps i>q in the prefactors gives 



T 



(584) _ 



(J! 



eq 
584 



-hv Q /T 



s, 



where 



y 



(114) 



(115) 



depends on the single parameter r c . We note that with- 



out the factor of £, Eq. (114 1 would be the standard 
decay rate formula with escape probability P csc — 1/rg 
(appropriate in the high optical depth limit). Thus £ can 
be thought of as a correction factor associated with the 
H I continuum opacity. 
Unfortunately, 



the integral in Eq. (115) is not 



amenable to direct calculation, since the integrand is ill- 
conditioned near y pa 0. We may obtain an alternative 
form by taking the integral from y = —Y to y = +Y (we 
will take the limit Y — > oo), plugging in Eq. (Ill I into 
the integrand, and then realizing that £(+V) — > due to 
our boundary condition, we find 



£ = £(-Y) + r c 



€(y) dy- 



(116) 



Formally, for r c > 0, we have £(— Y) —> for Y —> oo, 
but numerically one must keep the first term at small r c . 
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FIG. 9. The function £(t c ) (solid line) and our approximation 
(dashed line). 



Equation (111 I is stiff and can be solved by the back- 
ward Euler method; however even this is too slow to use 
in a "fast" recombination code. Therefore, we have con- 
structed a fitting function for B , valid to within 0.8% for 
all positive r c : 



£{T C ) « y/l + 7T 2 T C + 



7.74r c 
1 + 70r c 



(117) 



Note that this has the correct limiting behavior B —¥ \ 
for r c — > 0, and for positive r c has B > 1, as one would 
expect. The full function B (r c ) and the approximation of 
Eq. (117) are shown in Fig. 9] 

The use of Eqs. (114 1 and 1117]) arc sufficient to solve 
for the behavior of the 584 A line. 



an "on-the-spot" approximation: we assume that for each 
photon emitted in the 584 A line, there is a probability 
that the photon is re-absorbed of 



P, 



rcabs 



o -r] c (E [2lpo] ~E [23po] )fh 



(i 



(119) 



where the first factor is the probability that the photon 
can redshift from 584 A to 591 A without being destroyed 
by H I continuum opacity, and the second factor is the 
probability that, once it reaches the 591 A line, the pho- 
ton is indeed absorbed. Thus the 591 A line has two 
effects: first, it contributes to the net formation of the 
He I ground state in accordance with Eq. ( 118 ), and sec- 
ond it reduces the net formation from the 584 A line by 
a factor of 1 — P rca bs- Thus we write 



X (/Si 
1 reabs* 7 



1"591 \ 

e -(- E [23p"]-- E [llS])/ T 



(120) 



4- Hydrogen recombination corrections 



As a final step, we include in our rate equation the 
contribution to the ionization fraction due to changes in 
the hydrogen Saha equilibrium: 



r 



(H) _ 



dx 



HI 



dt 



(121) 



Saha 



which is obtained by two-sided finite differencing with 
Az = ±0.5. This contribution to x e is necessary for full 
accuracy at the very last stages of He recombination. 



3. The 591 A line 

The intercombination line He I] 2 3 Pf ->■ l x S at 591 
A can reach optical depths of order unity during helium 
recombination. Therefore it must be considered carefully 
However, because the damping wings are optically thin, 
it is only the line center that is of interest. Because there 
is negligible H I continuum opacity in the core during 
helium recombination, we use the Sobolev approximation 
for the 591 A line. If there were no 584 A line, then 
we would have blackbody radiation entering the line and 
could write the net decay rate as 



tt-hc 3 



(1 - e- T ^) - e~ iE ^^- E ^^ /T ^ , (118) 



where T591 is the optical depth of the 591 A line, and 
/jg-L = X[ 2 3po]/(3o:[ 1 i5]) is the equilibrium value of the 
photon occupation number in the line. Based on the 
ratio of Einstein coefficients [5*0] . we take T591 — 1.023 x 
10- 7 t 584 . 

However, the 591 A line can also absorb photons that 
redshifted out of the 584 A line. We treat this problem by 



5. Integration details 

We may now construct the overall electron fraction 
evolution, T = + J^ 584 ) + ^ 591 ) + F^) . The 

equation is stiff at early times, which is why we do not 
turn it on until the free electron fraction given by the 
post-Saha expansion differs from the Saha equilibrium 
value by more than 5 x 10 -4 . We then use the ODE 
integrator described in Appendix [C*) 

The treatment of helium is turned off once the abun- 
dance of He + ions is less than 10 per hydrogen atom, 
after which we switch to following hydrogen recombina- 
tion (we checked that hydrogen and helium recombina- 
tion never overlap and can be followed separateley) . 



D. Comparison to detailed calculations 

The ultimate test of the sweeping approximations 
made in this section is their comparison against more 
detailed computations of the He I recombination history. 
We compare against Switzer & Hirata [33J in Fig. 10 
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Helium recombination comparison 




1800 2000 2200 2400 2600 2800 3000 



z 

FIG. 10. A comparison of the HyRec He II— >I recombination 
history to that computed by the "full physics" code of Switzer 
& Hirata [331. The maximum deviation is 0.3%. 



bination histories for different cosmologies and construct 
a fitting function 6, 53]. Our point of view here is that 
the physics of primordial recombination is simple enough, 
and an exact calculation from first principles is now fast 
enough that there should be no reason to use fudge fac- 
tors and approximate correction functions. This is es- 
pecially relevant if one wishes to extend the standard 
recombination calculation by introducing "exotic" new 
physics. We would like to emphasize that the fast compu- 
tation presented here, using the EMLA method, is very 
well adapted for the computation of the recombination 
history, but that the standard MLA approach and fast 
interpolation methods may still be useful for the compu- 
tation of the recombination spectrum. 

We believe our code is accurate enough (aside from 
neglecting collisional transitions) and has a sufficiently 
small runtime to be incorporated in Monte Carlo Markov 
chains for upcoming CMB data analysis from the Planck 
mission. 



The maximum error is 0.3%, which is roughly equal to 
the stated theoretical uncertainty in the Switzer & Hirata 
calculation. 



VII. CONCLUSIONS 

Wc have presented a complete treatment of primordial 
hydrogen and helium recombination, including all the ef- 
fects that have been shown to be important so far. Our 
computation accounts for the multilevel character of hy- 
drogen and the non-equilibrium of angular momentum 
substates, radiative feedbacks, two-photon transitions, 
and frequency diffusion in Lya for hydrogen recombina- 
tion. For helium recombination, we account for HI con- 
tinuum opacity in the He I 2 1 P° — l x 5 line, decays in the 
2 3 P° — 1 X S intercombination line, and feedback between 
these lines. We have implemented all these effects in a 
single recombination code, HyRec, which can compute a 
recombination history in ~ 2 seconds on a standard lap- 
top for a given set of cosmological parameters. Provided 
collisional transitions can be neglected (which remains to 
be established), we estimate the errors of our computa- 
tion to be a few times 10~ 3 during helium recombination 
and a few times 10 -4 during hydrogen recombination, 
including both numerical errors and errors due to the as- 
sumptions and approximations made for physical effects. 
If collisional transitions are shown to have a significant 
effect on recombination, our code can be easily updated 
to account for them with very little loss of computational 
efficiency. 

It has been argued that corrections to the recombi- 
nation history due to radiative transfer effects are rela- 
tively independent of cosmology [52] , and that one could 
therefore compute them once and use the resulting cor- 
rection function to account for them for any given cos- 
mology. Alternatively, one might run a grid of recom- 
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Appendix A: Proof of some relations involving 
effective rates 

1. Preliminaries 

Here we use the same notation as in Paper I. Capital 
indices K, L refer to "interior" excited states, and lower- 
case indices i,j refer to "interface" excited states. We 
define the rate matrix M whith coefficients: 



M 



KL 



?KLt K 



(1 



>KL 



(Al) 



In Paper I, we have shown that the populations of the 
interface states are given by: 



X K = £ (M _1 ) 



LK 



(A2) 



We also showed that the probabilities P % K and P^ are 
given by 



and 



P* K = J2(M- 1 ) KL RL,i, 

L 

p* K = j2(M-i) KL h. 



(A3) 



(A4) 
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In Appendix C of Paper I, we showed that M satisfies 
the following detailed balance relation: 



Qk {m- 1 ) kl = q l (m- 1 ; 



LK ■ 



(A5) 



where Qk = 9k^~ Ek ^ Ti is the contribution of individual 
states to the partition function and gx is the degeneracy 
of the state K. 



2. Rate of change of the free electron fraction 



In this section we derive Eq. (31), which was not de- 
rived in Paper I. In the standard MLA formulation, the 
rate of change of the free electron fraction can be written 
as: 

x e = -y~] [n H x 2 e a K - X K /3 K ] 

K 

- ^2 [nnxtati - . (A6) 

This formula is never used in standard MLA codes, as 
it requires a summation over a large number of nearly 
cancelling terms, and MLA codes use instead x e — —xi s 
to compute the rate of change of the free electron frac- 
tion. Eq. ( |A6[ ) remains however formally correct. Using 
Eqs. (A2) and (A4|, we rewrite: 



LK 



K 



K.L 



H,L 



= Y P l n n x 2 e a L + ^2 x i R i 

L L t 

= nnx 2 e a L - n H x 2 ^ E a L P l L 

L i L 



(A7) 



where in the last equality we have used the complemen- 
tarity relation P K +P K = 1- Inserting this result into 
Eq. ( A6 1 , and using the definitions of the effective recom- 



bination coefficients and photoionization rates Eqs. ( 23 ) 



and (24), we immediately recover Eq. (31) 



3. Proof of Eq. (41) 



Consider Eq. £6|) with T K m T K (T r ). The formal 
solution for the Pjf is given by: 

P}f^J2( M ^)KL^- (A8) 

L 

Therefore one may rewrite Eq. ( [34| , for i = 2s, 2p: 

Hi,ls = Ri,ls + E ^i,K(T T )RK,ls, (A9) 
K 



where we have defined 

X itK (T T ) = J2Ri,L (M- 1 

L 

Qk 



LK 



1 KL 



5^ e -^/T rp ^ (Tr)j 



(A10) 



where in the second line we used Eq. ( A5 1 , in the third 



line we used the detailed balance relation verified by Ri l 
and R.L.ii an d in the last line we used the for ma l solution 
for P K , Eq. (A3). We therefore obtain Eq. (41). 



4. Expression of Xk in terms of Xi,x e and effective 

rates 

Taking T m = T r and using the detailed balance relation 
j e jiH«L = QlPl, we rewrite Eq. (A2| as follows: 



Qi 



LK 



Rl, 



L 

+ E^E( M ' 1 )kl^ (AH) 



where in the last equality we have used Eq. ( A5 ) . Using 
the formal solutions for the probabilities Eqs. (A3), (A4), 



we see that we recover Eq. ( 42 ) 



Appendix B: Extrapolation of the effective rates to 



We have tabulated the effective rates 
A- 2s {T m ,T t ),A2p{T la ,T T ) and TZ2s,2 P (Tr), including 
all excited states up to the principal quantum num- 
ber rt max , for several values of n max up to 600, over 
the temperature range 0.004 eV < T r < 0.4 eV, 
0.1 < T m /T r < 1. This range of temperatures corre- 
sponds to 20 < z < 1650 for a wide range of cosmologies. 
For every pair (T m , T r ), we have fitted the effective rates 
by the following functional form: 



Ai{T ni: T t \ 7i max ) — l A.j / (T rn: T F ] oo) I 1 



(Bl) 

and similarly for lZ2s,2p, where k and 7 depend on T m 
and T T as well as on the coefficient being fitted. This 
allows us to extrapolate the effective rates to n max — > 00. 
Of course, this is only a formal extrapolation, as for 
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n larger than a few thousands, the excited states of 
hydrogen are no more well defined (see Ref. [54] for 
a discussion). The extrapolated rates are still more 
accurate than those computed with a a finite number 
of states. The residuals of the fit have a maximum 
relative amplitude of 5 x 10 -4 over the whole range 
of temperature considered, for 200 < n max < 600, 
and more than an order of magnitude smaller on the 
restricted range T t > 0.04 eV, T m /T r > 0.8 which 
corresponds to z > 200 (note that neglecting the 
overlap of the high-lying Lyman lines leads to errors 
in the effective rates of similar amplitude [T5]). For 
reference, the maximum relative difference between the 
effective rates computed with n max = 600 and their 
extrapolation at n max = oo is 0.05 over the whole range 
of temperature considered, and 0.002 over the restricted 
range corresponding to z > 200. We checked that 
our method recovers the correct case-B recombination 
coefficient a B (T m ) = Y*=2b,2 P A(^m, T r = 0;oo). Our 
extrapolated a b agrees with the fit of Ref. [55] to better 
that 0.2 % for T m > 40 K, which is the accuracy claimed 
by the authors of Ref. [55] . 



Appendix C: Numerical ODE integrator 

For the sake of computational efficiency, we use a sec- 
ond order ODE integrator that uses derivatives computed 
at previous timesteps. This allows us to evaluate deriva- 
tives only once at each timestep. Explicitly, to numer- 
ically solve the equation y'(x) = f(x,y), we use evenly 
spaced steps Ax, and obtain the solution at the (i + l)th 
step as follows: 



Appendix D: Post-Saha expansion at early phases of 
hydrogen and He II^I recombinations 

As explained in Appendix D of Paper I, the ODE 
describing hydrogen recombination is stiff at z > 1500 
and so is the ODE describing He II— > I recombination at 
z > 2800. We therefore use an expansion around the 
Saha equilibrium solution: 



x° + 




(Dl) 



where xf is the Saha equilibrium value of the free electron 
fraction. 

1. Hydrogen recombination 

The Saha equilibrium value of the free electron fraction 
is the solution of the following equation: 







n S\2 



--E//T 



(D2) 



where g e was given in Eq. ([3]) and T = T m = T r at 
early times. The numerator in Eq. (Dl I can be obtained 



analytically by differentiating Eq. ( D2 ) 



d(x s e ) 
dt 



#(f -§)(* 
2zS + 



S\2 



(D3) 



For the denominator in Eq. (Dl), we numerically differ- 
entiate the derivative x e obtained when accounting for 
two-photon processes and diffusion, using a two-sided fi- 
nite difference with Ax e — ±0.01(1 — a;f). 



He II^I recombination 



Vi+i = Vi + &Vi, 
A Vi = Ax [1.25J/J 



0.25yU 



(CI) 



where y[ — f(xi,yi) is stored at each timestep for later 
use. For the case of interest, we have x = In a, y = x e 
and / = x e /H . 



In that case the free electron fraction in Saha equilib- 
rium is given by x^ — 1 + q, where q can be obtained 
from Eq. (97). As in the hydrogen case, differentiation 



of Eq. (97) gives us an analytic expression for the nu- 
merator in Eq. (Dl). We numerically differentiate the 
derivative x P given by Eq.(99) using a two-sided finite 

in- 



difference with Ax e — ±0.01(1 + /n e 
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